Event-Driven Dynamics of Rigid Bodies Interacting via Discretized Potentials 
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A framework for performing event-driven, adaptive time step simulations of systems of rigid bod- 
ies interacting under stepped or terraced potentials in which the potential energy is only allowed 
to have discrete values is outlined. The scheme is based on a discretization of an underlying con- 
tinuous potential that effectively determines the times at which interaction energies change. As 
in most event-driven approaches, the method consists of specifying a means of computing the free 
motion, evaluating the times at which interactions occur, and determining the consequences of in- 
teractions on subsequent motion for the terraced-potential. The latter two aspects are shown to 
be simply expressible in terms of the underlying smooth potential. Within this context, algorithms 
for computing the times of interaction events and carrying out efficient event-driven simulations are 
discussed. The method is illustrated on system composed of rigid rods in which the constituents 
interact via a terraced potential that depends on the relative orientations of the rods. 
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I. INTRODUCTION 

Simulating systems with discontinuous interactions of- 
fers a number of advantages over standard molecular dy- 
namics simulations (SMD) in which the solution of a sys- 
tem of ordinary differential equations is solved numeri- 
cally by iterating a map which approximates the short- 
time dynamics[l|. The advantage of simulating discontin- 
uous systems using event-driven algorithms Q (discontin- 
uous molecular dynamics, or DMD) over SMD is particu- 
larly apparent for systems of low density, where the short- 
time mapping of the dynamics in SMD is applied to freely 
evolve the majority of the system. In spite of the rar- 
ity of interactions at low density, the inherent time step 
in SMD simulations cannot be increased beyond some 
small threshold without loss of stability and accuracy of 
trajectories, since there is nearly always a small fraction 
of particles interacting at any given time. Several adap- 
tive integration methods exist, typically based on either 
separating rapidly and slowly-varying components of the 
potential in a multiple time step approach[3[, or on time 
re-parametrization of the Hamiltonian equations to a new 
system that is integrated with a fixed step sizeQ- Both 
methods have inherent drawbacks making them unsuit- 
able for arbitrary potentials. On the other hand, in the 
event-driven approach where there is no inherent time 
step, each non-interacting particle is not propagated for- 
ward in time until it interacts with another particle in 
the system. 

Event driven simulations also have some rather serious 
drawbacks. For example, simulations of flexible molecu- 
lar systems are plagued with processing often irrelevant 
intra-molecular events on very rapid time scales, wasting 
a great deal of computation. When physically reasonable, 
much can be gained by treating the molecules as rigid. 
Although the general framework for performing event- 
driven simulations of rigid or constrained systems has 
been recently worked out[5|, numerical methods must be 
used to find interaction times, potentially leading to gross 
inefficiencies. Another clear drawback that discourages 



the use of the DMD approach is that little is known about 
how to design site-based, stepped interaction potentials 
between such rigid bodies, and much work must be done 
to tune interaction parameters, such as well-depths and 
interaction distances at which discontinuities occur. The 
design of a distance-based, discontinuous potential for 
long-ranged electrostatic interactions seems particularly 
problematic. 

In this article, possible solutions for both these prob- 
lems in DMD simulations are presented. An algorithm 
based upon an adaptive grid search is presented as an al- 
ternative to the uniform grid approach of Ref . [H to make 
the search for interaction times as efficient as possible. 
We also demonstrate here how the issue of designing de- 
tailed discontinuous potentials can be side-stepped alto- 
gether by using a mapping of an underlying continuous 
pair potential onto discrete potential energy values. The 
potential energy then consists of a set of allowed energy 
terraces, each mapped onto by many different positions 
and orientations of the system. The discretization of the 
potential energy on the level of the pair potential implies 
that the evolution consists of free propagation of the sys- 
tem punctuated by impulses at discrete times when the 
underlying continuous interaction potential for a pair of 
particles hits a critical value. Because the scheme uses 
ordinary continuous interaction potentials as its basis, no 
substantial effort is required to parameterize the Hamil- 
tonian, and the experience of many years of work in the 
modeling of systems can be exploited. The method is 
applicable to any type of pair-interaction potential, and 
can be used with potentials written for rigid systems that 
depend on the center of mass positions as well as the rel- 
ative orientation of two interacting bodies [1,0]. 

The paper is organized as follows: Sec. HT1 reviews the 
elements of rigid body mechanics required to formulate 
the method. In Sec. IIIH a general description to numeri- 
cally find interaction times is presented. A derivation of 
the consequences of the action of the impulsive forces and 
torques on the system is given in Sec. HVl The scheme 
is applied to a simple model system with orientationally- 
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dependent interactions in Sec. [V] Final comments are 
given in Sec. IVI1 



II. RIGID SYSTEMS INTERACTING VIA 
STEPPED-POTENTIALS 

The systems considered here consist of N rigid bodies, 
each of mass m and moment of inertia tensor lj. Associ- 
ated with each object i are a center of mass position and 
velocity and Vj, an attitude or orientation matrix A i; 
and an angular velocity vector u>i. The attitude matrix 
Aj transforms a vector a with respect to a fixed inertial 
lab reference frame to its representation a = A.^a in the 
principal- axis frame of body i. Aj is in fact defined such 
that the moment of inertia tensor is diagonal in the prin- 
cipal axis frame: 1; = AjljAj = diag(iji , U2, Uz), where 
lit, Ii2 and Ji3 are the (possibly distinct) principal mo- 
ments of inertia of body i. Although there are several 
ways to parametrize the attitude matrix A, such as Eu- 
ler angles, unit quaternions, or angle axis 8 coordinates, 
three generalized coordinates are always required to spec- 
ify the orientation of each three-dimensional rigid body, 
denoted here by #j = (flu, ^3). The time derivative 
of Di can be related to by noting that Wj is related to 
the time derivative of the attitude matrix viajlcj 

^2 £bac(Ui) a = ^ (Ai)ab(K)ac, (1) 

a—x,y,z a—x : y : z 

where £(, oc is the Levi-Civita symbol, [s] From this re- 
lation, one can easily derive that = N|i?j, where 
(Ni)ah = \e bcd {k i ) ec d{k i ) ed /d{'d i ) a . 

If the system is governed by a smooth potential U, 
then the equations of motion imply that 

where pj = mv,; is the (linear) momentum of body i, 
hi = \iU>i is the angular momentum of that body with 
respect to its center of mass, Pi is the force on the center 
of mass of body i, and Tj is the torque on body i. 

In many cases, the formal expression for the torque in 
Eq. ^ can be written in compact form which does not 
depend on the choice of parametrization of the attitude 
matrix Aj. For example, if the potential can be written 
in terms of site-site interactions, the center of mass force 
Fj = Fi 7 is a sum of the forces Fj 7 acting on the 
sites 7, and the torque can be written as 

n = ^2 ( r n - r i) x F *7' ( 3 ) 

7 

where r^ 7 is the position of site 7 on body i. 

A second example where one can write simple expres- 
sions for Fj and r, is if the potential U is a sum of pair 
potentials Uij between bodies i and j, each of which is 
rotationally invariant but depends on inner products of 



the relative position vector = r-j — rj and a set of ori- 
entationally dependent vectors s" and (where a and 
j3 are integers indicating different vectors for body i and 
j, respectively). Then the force and torque on body i 
due to interactions with body j through the interaction 
potential [/^-(ry, {sf}, {s^}) can be written as 

(4) 
(5) 

where we have used Eq. ((Sj) and the fact that ebde(^i)ab — 
(Aj)cd d(Ai) ce /d('di) a . Using the expressions of the forces 
in Eq. Q , one finds for each interacting pair Fy + Fji = 
0, which implies conservation of total linear momentum 
Pi- Furthermore, using the torques in Eq . ([5]) . it is 
straightforward to verify that 

r s ; X Fy + Tj X Fji + T tJ + T 3i = 0, (6) 

which, in turn, implies that the total angular momentum 
X^(Li + fj x pj) is conserved by the dynamics. fill] 

So far we have considered U to be a continuous po- 
tential, which is required for the derivatives in the above 
equations to exist. Now we consider a stepped interaction 
potential between a pair i, j of bodies bases on a con- 
tinuous potential Uij , such that the interaction potential 
between bodies i and j is of the form 

K 

Vij = V mm + - U k )AV k> (7) 

k=l 

where Q(x) is the Heaviside function, and Uk are a dis- 
crete sets of values of the smooth potential at which the 
system gains or loses potential energy AT4. Note that 
for a system in which the underlying continuous interac- 
tion energy Uij between bodies i and j lies anywhere in 
the range Uk < Uj < Uk+i, the form of Eq. (J7J) assigns a 

constant potential energy value of V m in + J2k'=i ^^fc' to 
this interaction. Figure [1] contains an illustration of the 
result of this procedure for the potential energy function 
given in Eq. I]22p. Because of the shape of the potential 
energy landscape, we call a potential of the form Q a 
terraced potential. 

For a system in which all interactions are of terraced 
form, the procedure to evaluate dynamical properties is 
the same as that in simulations of hard sphere systems[5|. 
While the system has constant potential energy between 
interaction events, the motion of the system is free and 
the trajectory of each constituent is independent of all 
others. The free propagation of the system between 
events determines the evolution of the spatial coordinates 
of the molecules in the system, and the interaction times 
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FIG. 1: Illustration of the potential energy discretization pro- 
cedure given in Eq. (0 applied to the potential in Eq. (|22p 
of Sec. [V] with the e = 5/2, a = 1 and A = 1. The top 
panel shows the continuous potential U and the bottom panel 
the terraced variant V if Vmin = — 2, A 14 = AV — 1 and 
Uk — Vmin + (k — \)AV. On the axes, r = \n — and 

C = Si • Sj . 



at which the momenta of the system change discontin- 
uously is determined by identifying the times at which 
the underlying continuous pair potential energy Uij hits 
an energy terrace Uk where the potential energy changes 
discontinuously. Even though this interaction time is de- 
termined by free motion of the system, in general, it must 
be found numerically due to the mathematical complex- 
ity of the interaction condition. The next section ad- 
dresses this problem. The final ingredient required to 
perform event-driven simulations of systems consists of 
specifying the interaction rules of how the momenta of 
the constituents in the molecular system are altered at 
the interaction times and is worked out in Sec. IIVI 



III. FINDING INTERACTION EVENT TIMES 

For systems with terraced potentials as in Eq. ([7]), the 
particles evolve freely until a configuration of the system 
is reached where the instantaneous value of the underly- 
ing continuous potential fZy equals Uk, for some pair of 
particles i and j. The first time at which such an event 
occurs can be solved by finding the smallest positive zero 
(root) of the indicator equation 

f l3 {t) = U l3 (v t] {t), K(i)}, {sj(t)» - U k , (8) 

where the time dependence of the indicator function 
is determined by the free trajectories of bodies i and 
j. Although the time dependence of the relative vec- 
tor Vij{t) = Tij(0) + (v, —Vj)t appearing in the indicator 
function is simple in the case of free motion, the time- 
dependence of the orientational vectors sf (t) and Sj(t) 
depends on the evolution of the attitude matrices A^(i) 
and Aj (t) , since these vectors may be expressed as 

sf(t) = At(t)if (9) 
sf(i)=A+(i)sf, (10) 

where sf and are time-independent vectors in the 
body frames of particles i and j, respectively. 

The form of the time dependence of the attitude matri- 
ces Ai(t) depends on the way in which mass is distributed 
in the rigid bodies. Nonetheless, it is possible to write 
down analytical expressions for the time-dependence of 
Aj even when the mass of a body is not symmetrically 
distributed. The general solution of Aj(t) can always be 
written in the form[Io| 

A i {t) = P i (t)-A i (0), (11) 

where the matrix Pj (t) propagates the orientation matrix 
from an initial time to time t. The precise forms of Pj(t) 
for different kinds of rotor can be found in Ref. [Tol . 

Even with analytical expressions for the time- 
dependence of the underlying potential Uij(t), the ear- 
liest interaction time must be found numerically, as was 
also the case for event-driven dynamics for rigid systems 
interacting via site-site potentials Q. However, one ma- 
jor benefit of utilizing an energy terracing approach is 
that there is only a single one-dimensional root search to 
be conducted for each pair of interacting bodies, in con- 
trast to a site-based energy approach in which all pairs 
of sites between pairs of molecules must be examined for 
the earliest interaction time. 

In Ref. a simple approach to find the earliest in- 
teraction time was given in which screening methods are 
used to identify a minimum and maximum time between 
which a root of the indicator function could lie. This 
interval is then sub-divided into equally-sized smaller in- 
tervals of fixed size At and used to bracket sign changes 
of the indicator function. Unfortunately, for translating 
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and rotating rigid bodies, the indicator function is oscil- 
latory, making the detection of so-called grazing interac- 
tions troublesome unless a very small grid interval At is 
used. Since the indicator function has a local extrcmum 
in a grazing interaction, a good strategy to find this kind 
of event is to determine the minimum or maximum of the 
indicator function in cases in which the indicator func- 
tion fij(t) itself does not change sign, but its derivative 
fij(t) does. However, in the vicinity of a grazing inter- 
action, the values of the indicator function are typically 
small, and hence only extrema where the indicator func- 
tion at one of the grid points lies below some thresh- 
old value need to be investigated. Efficient numerical 
routines exist to find local extrema of a one-dimensional 
function (l^. [l3l| . and thereby detect grazing interactions. 
Once an interval with a change in sign of the indicator 
function has been identined,standard techniques can be 
invoked to find the root [13, El • 

Although the scheme outlined above is fairly robust 
and can detect many millions of events without missing 
roots, its efficiency is strongly dependent on the choice 
of basic time interval At. In many cases, the size of this 
interval is determined not by the typical rate of change 
of the indicator function but by certain rare scenarios 
where rapid changes in the indicator function and its 
derivative and multiple local extrema occur. Unneces- 
sarily small time intervals can be avoided by using the 
following adaptive bracketing scheme based on cubic in- 
terpolation to estimate the indicator function and the 
number of extrema in a given interval. 

The basic idea is to use information from previous 
grid points as much as possible. Note that when eval- 
uating the indicator function / at a certain point t, it 
is very easy to also compute its time derivative. For a 
given pair of molecules i and j and assuming a terraced- 
interaction potential based on the continuous pair poten- 
tial £/y(ry , {s.i} , {Sj}), the time derivative of the indica- 
tor function /y (t) is given by 



/«(*) 



(12) 



where Vy = Vj — Vj is the relative velocity vector for the 
centers of mass, Fy is the force on i due to j, and Ty and 
Tji are the torques on i and j, as given by the smooth 
interaction in Eqs. Q and ((5J. 

At the start of the search for the first root of the func- 
tion fij, its value and derivative are computed. A linear 
extrapolation fa(t) = fi + fit can be used to get an es- 
timate t' for the root of /y by solving fe(t') — 0. The 
next grid point is then taken to be t% = t' + St, where 
a small St is added to enhance the probability of a sign 
change of fij in [0,f2]- At t — t2, the function and its 
derivative are evaluated. If there is a sign change in /y , a 
root has been bracketed and a numerical root search us- 
ing standard techniques is used [12, HI] . Otherwise, given 
the value of the indicator function and its derivative at 
two times, a unique cubic interpolation of the indicator 
function can be constructed. For example, consider the 
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FIG. 2: Example of the adaptive root search with cubic in- 
terpolation. 



indicator function f(t) satisfying /(0) = fx, f(At) = f-2 
with time derivatives /(0) = /i and /(At) = f^. Us- 
ing these values, the cubic approximation f c {t) for the 
indicator function over the interval is 



fc(t) = fi + fit + at 2 + l3t 3 



where 



= 



3(/ 2 - /i) - (2 A + / 2 )At 
At 2 

2(/i - h) + (A + A)At 
At 3 



(13) 

(14) 
(15) 



The number of extrema in the interval [0, At] is then esti- 
mated by finding the number of real roots of the equation 
fc(t) — in the interval. Furthermore, the values of the 
extrema are easily obtained using (|13[) . 

The adaptive procedure using the cubic approximants 
is most-easily explained by considering the example in 
Fig. [21 Consider beginning the process of looking for 
a root after an interaction event. In this scenario, at 
point (1) in the figure, only the value of the indicator 
function and its derivative are known at t = 0. The next 
bracketing point is chosen by solving fi(t 2 ) = 0. If the 
time t' 2 is very large or infinite, the next bracketing time 
t<i is chosen to be some default (and regular) value, At. 
On the other hand if t' 2 < At, the second bracketing time 
ti = t' 2 +8t is chosen. At point (2), the indicator function 
and its derivative are evaluated and the cubic and linear 
approximations to the function evaluated. The linear 
and cubic equations are then examined for roots. In the 
case shown in Fig. [51 the cubic approximation has no 
roots whereas the linear interpolant does, at t' 3 . Since 
t' 3 < t<x + At, the next bracketing point is selected to 



be U = t'o 



St. But the cubic interpolant has a local 



minimum in the interval [£2,^3], so a minimum search 
algorithm is used to find the actual minimum of /y at 
(3). Since /y is found to be positive at (3), a potential 
grazing interaction has been ruled out. The root of the 
cubic t\ is then used to place the next bracketing point 
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at ti = t± + St. The cubic interpolant using points £3 and 
ti reveals that there is a local maximum in this interval, 
and that there is a sign change in the indicator function 
in the interval since f(ti) < 0. The local maximum is 
investigated and found to be positive at time t§, so the 
final bracketing interval where a root is found is taken to 
be [£5, £4], leading to the root at point (6). 



IV. INTERACTION RULES 



furthermore gives 



dUi. 



x s; 



dUi, 



+ (S\ - S f ) 



dUi. 



i) 



(s?xr«)=0. (20) 



If the total interaction potential for the system is a 
pairwisc-additivc combination of terraced two-body in- 
teraction potentials of the form ([7]), then forces and 
torques only act instantaneously on a pair of molecules i 
and j at a time determined by when [Zy is equal to one 
of the reference values Uk- The forces and torques are 
therefore impulsive, and lead to discontinuous changes 
in the linear and angular momenta of the pair. The ef- 
fect of the impulses must be consistent with conservation 
conditions that arise from the symmetries of the over- 
all Hamiltonian, such as the conservation of energy, and 
conservation of total linear and total angular momenta. 

Since the momenta and angular momenta change dis- 
continuously and are not well-defined at an interaction 
time t c , a more practical starting point for deriving the 
interaction rules is to consider the effect of impulses ap- 
plied to the overall change in the momentum and angular 
momentum of interacting bodies i and j over a small time 
interval [t c — 6, t c + S] : 



(16) 
(17) 



rtc+S 

p< = p 4 + / Fi(t)dt 

Jt c -S 
rt c +S 

14 = U+ / n(t)dt, 

JUS 



with analogous expressions for molecule j. In Eqs. |T| 
and l|17p. the primed and unprimed vectors represent 
post and pre-interaction values of the respective quan- 
tities. For small enough 5, the probability of a parti- 
cle other than j interacting with i becomes zero. From 
Hamilton's equations for the discontinuous system, the 
impulsive forces and torques are then given by 



F<(*) 



(t c )6(t - t c ) 



SlnjiQSit-Q, 



(18) 
(19) 



where Py(t c ) and 7y(i c ) are the forces and torques on 
body i at the interaction time t c for the continuous sys- 
tem given in Eqs. (|3]) and ([5]), and S- and S? are unknown 
scalars. Note that the directions of the impulsive forces 
and torques due to the pair interaction are along the 
directions of those quantities for a continuous system in- 
teracting by the same potential at the interaction time. 
For the interaction pair i-j, conservation of linear mo- 
mentum immediately implies that Sj — S? — . The 
requirement of conservation of total angular momentum 



Since this condition must be satisfied for arbitrary vec- 
tors r ih sf, s^, we conclude that = Sj = S f = 5.Q 
The unknown scalar S is now found from conservation of 
energy, by solving the quadratic equation 



l'\ I\ , T, • I, ' • T, 
+ % + 

m 2 

+ (v y - ■ Fi + n ■ U>i 



S' 2 



•w j )5 + Ay = 0, (21) 



where all quantities are evaluated at interaction time t c 
and AV is the change in potential energy. Note that the 
term proportional to S can also be written as —Uij (t c ), as 
follows from Eq. (p~2|) . The physical solution corresponds 
to the positive (negative) root branch if Uij(t c ) < 
(Uij(t c ) > 0), provided real roots exist. If this latter 
condition is not met, there is not enough kinetic energy 
to overcome the discontinuous barrier, and the system ex- 
periences a reflection with no change in potential energy, 
i.e., S is the non-zero solution of Eq. (|2ip with Ay = 0. 



V. MODEL SYSTEM 

As an illustration of the method, consider a system 
composed of rods in which the continuous interaction po- 
tential between a given pair i-j is of a modified Lennard- 
Jones form 



U i:i = 4e 



12 



with 



= |r«|+C 



- (s 4 • s^ 5 



(22) 



(23) 



Here, Sj is a vector pointing along the long axis of the 
rod and r e defines an orientation dependent effective dis- 
tance of interaction. Note that the form of this potential 
makes it energetically favorable to align adjacent rods 
orthogonal to one another. 

By construction, the potential is invariant under ro- 
tations and translations and therefore the dynamics con- 
serves the total linear and angular momentum in addition 
to the total energy of the system. Using Eqs. (31) and 
([5l), the force and torque on body i due to interactions 
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with j is 



24e 



(24) 



Note that the torque on i is orthogonal to the vector s,-, 
so that the component of the angular momentum along 
the axis of the molecule is constant. Thus, even though 
the rotational dynamics of the rigid rods is that of a sym- 
metric top, the angular motion around the long axis is 
decoupled from the other directions and can be neglected. 

From this continuous system one can construct a sys- 
tem with a terraced potential. In the current study, the 
simplest form of the terracing procedure will be used, 
where the discontinuities in the potential are evenly- 
spaced values of the interaction energy, i.e. in Eq. ([7]) 
all AV4 are taken to be equal and the values of Uk are 
taken to be Uk — (k — |)AV. in practice, it is desirable 
to map a value of w to zero, so that pairs that are 
far apart do not contribute potential energy to the sys- 
tem. This is readily accomplished by adjusting the value 
of Vmin in Eq. (|7|) . It is also important that the critical 
values Uk not coincide with minima of the smooth pair 
potential, because a terrace of width zero would result, 
leading to numerical instabilities. 

It is interesting to consider how the graininess of the 
energy terraces influences dynamical phenomena, par- 
ticular in dense systems where orientational relaxation 
times become quite long. To examine the sensitivity 
of dynamical correlations to the form of the stepped- 
potential, event-driven simulations of a system of N = 
512 rods at a reduced number density p* — pa 3 = 0.5 and 
reduced temperature T* = kT/e = 1.4 were carried out 
in a periodic, cubic box and contrasted with SMD simula- 
tions using the continuous underlying potential, Eq. (|22j) . 
The parameters chosen for the system and pair potential 
were e = 3.0 kJ/mol, a = 3.035 A, ( = 1 A, m = 18 
g/mol, and h = h = 1.154 g/mol-A 2 . 



A. Implementation of the event-driven simulation 

As was observed in the case of event-driven simula- 
tions based on site-site potentials Q, the efficiency of the 
event-driven approach can be greatly enhanced by im- 
plementing several simple and now fairly standard tech- 
niques, such as cell divisions, tree data structures to 
manage events, and local molecular clocks^. In Ref. [f| 
two additional techniques for improving the performance 
of event-driven simulations utilizing numerical methods 
of finding event times were also presented: the use of 
screening methods to identify initial and final bracket 
times, and the truncation of potentially non-useful event 
searches through the scheduling of a virtual interaction. 




FIG. 3: The radial distribution function for the center of mass 
of the rods. 



To identify which pairs of particles could have interac- 
tions under the terracing potential for a given choice for 
the set of discontinuities, the maximum distance r m at 
which U is larger than AV/2 was found by solving the 



equation L7y = —AV/2 with r e — 



C/2, giving 



i[l-Vl-AV7(2e)] 



(25) 



The system is then partitioned using this maximum in- 
teraction distance to determine the cell size so that only 
particles in the same or adjacent cells interact. For each 
pair of possibly interacting particles i and j, the times at 
which the distance between center of masses reaches r m 
is solved exactly using the linear free motion of the pair 
to determine minimum and maximum bracketing times. 

To implement the virtual interaction event, a root 
search procedure for a given pair of rods was only con- 
ducted up to a maximum of 4 femtoseconds before the 
root search was truncated, and if no root was found, the 
value of the indicator function and its derivative at the 
end point of the interval were stored in the tree to resume 
the search for the pair if neither particle has an event 
before the point at which the search was halted. A max- 
imum step size of At = 4 femtoseconds was also utilized 
in the adaptive search algorithm outlined in Sec. IIII1 

The use of a uniform maximum step size generally 
leads to sub-optimal tree structures in the priority queue, 
where long linear branches corresponding to scheduled 
simultaneous virtual interaction events lead to longer in- 
sertion times of new events into the queue. Two coun- 
termeasures were employed to improve the management 
of the event queue: First, the maximum step size At 
was given a small random adjustment to avoid the co- 
incidence of any two virtual interactions. The effect of 
this adjustment is to form more balanced binary trees 
free of long linear branches. Second, a bounded prior- 
ity queue [15j was used in which arrays of linear lists of 
events in a given time interval are combined with an im- 
plicit heap binary tree that performs a fine sort of the list 
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containing the current time interval. The use of the lin- 
ear lists enables rapid insertion of future events into the 
priority queue, and typically reduced the overall compu- 
tational time by around 15% for the system sizes investi- 
gated. It is expected that the advantage of the bounded 
priority queue over standard binary tree data structures 
will increase with increasing system size, since operations 
on the queue take 0(1) time per event, as opposed to 
0(log N) time. 

To keep the tree structure from growing too large in 
the course of a simulation, standard scheduling schemes 
require canceling each event which becomes invalid due 
to the occurrence of an earlier interaction involving one 
of its participants. Tracking down these events in the tree 
requires a fair amount of bookkeeping Q. Directly can- 
celing invalidated events could be done with the bounded 
priority queue as well, but there is no real need to do so 
since the tree component of the event queue is small al- 
ready. Not removing invalid events from the queues was 
found to lead to another 15% reduction in computational 
time. The problem of growing linear lists was handled 
by a (relatively fast and infrequent) cleanup of invalid 
events from the lists when computer memory threatens 
to be depleted. 



B. Comparison between standard molecular 
dynamics and event-driven simulations 

To assess the relative merits of the DMD approach as 
opposed to a standard rigid body molecular dynamics 
approach, SMD simulation were carried out based on the 
pair potential in Eq. (|22p . The equations of motion were 
integrated using a symplectic integration scheme that uti- 
lizes the exact free motion for both rotational and trans- 
lational degrees of freedom The time step used in 
the SMD simulations was 1.8 femtoseconds, resulting in 
relative fluctuations of the total energy to the potential 
energy of about 1%. To improve the efficiency of the 
SMD simulations, Verlet lists were used.[l| In addition, 
the potential was smoothly interpolated to zero starting 
from a cut-off distance of ri — 2.5a and reaching zero 
at a distance of r u = 3<r, by taking as the interaction 
potential 



smooth 



9(nj)Uij(r. 



l 3 V 1J i * ' 3 / 



where 



1 r <v\ 

g(r) = ^ ^ — {ru _ ri) * ' n<r <r u 
r > r u 



(26) 



(27) 



The effect of discretization in the energy-terrace model 
on static structural and dynamical correlations was ex- 
amined by varying the energy difference between discon- 
tinuities from a fraction of kT, the natural energy unit of 
the system, to several multiples of kT. Not surprisingly, 
the radial distribution function for the center of mass of 




0.5 0.75 
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FIG. 4: The normalized center of mass velocity autocorrela- 
tion function for different levels of energy terracing. 



the rods, shown in Fig. [3l shows the same qualitative 
behavior as that obtained from the continuous potential 
system, with clear discretization effects visible when large 
energy gaps (i.e. AV = kT) are used. As observed in 
other systems with discontinuous potentials @, the peaks 
and troughs observed in the radial distribution function 
tend to be exaggerated, but generally in a manner that 
results in the correct integrated number of neighbors for 
each radial shell. The radial distribution functions from 
the SMD and DMD simulations are in excellent agree- 
ment when AT/ is kT/2 or less. 

Dynamical correlations, on the other hand, tend to 
be more sensitive to the level of discretization in the 
stepped-potential. Considering the normalized autocor- 
relation function of the velocity of the center of mass 
(VACF), shown in Fig. gj it is clear that the VACF de- 
cays too quickly for large steps. This behavior can be un- 
derstood by noting that the average number of neighbors 
in the first solvation shell is too large when AT^ = kT, 
as can be seen from the radial distribution functions in 
Fig. [3] The increase in the average number of neighbors 
around a given particle leads to an enhanced number of 
interactions at short times that give an exaggerated kick 
in random directions to the center of mass velocities and 
leads to a rapid loss of correlation. However as the energy 
steps are reduced, there is a convergence of the VACF to- 
wards the result of the SMD simulations since the local 
structure as well as the magnitude of the impulses are 
modified. Note that when AV^ = kT/2, the life-time of 
velocity correlations is roughly correct. 

The same trends are observed in the orientational cor- 
relation functions. As is clear from Fig. [5], the normalized 
autocorrelation function of the long-axis of the rods de- 
cays too slowly when the energy discontinuities are too 
large. Once again, the long lifetime of orientational cor- 
relations for the system with large discontinuities is in- 
dicative of pairs of molecules being too tightly bound in 
an orientationally dependent, locally- preferred configu- 
ration. However, the DMD results rapidly approach the 



8 




FIG. 5: The normalized autocorrelation function of the long 
axis vector for different levels of energy terracing. 



SMD result as AV is reduced to kT/2 and exhibit the 
correct degree of anti-correlation for times around 0.15 
ps. 

For dense systems, the event-driven simulations are not 
expected to be much more efficient than SMD simulations 
that take of advantage of stable and accurate symplectic 
integration schemes. The cost of the DMD simulation 
scales linearly with the number of events processed per 
picosecond of time. Interestingly, for AV in the range of 
kT to 2 kT, the number of interactions per picosecond 
is approximately constant and equal to 3.5 x 10 4 coll/ps. 
At this rate of events, the DMD simulations are roughly 
1.7 times more efficient than the SMD simulations. How- 
ever for AV = kT/2, where the dynamics compares well 
with the dynamics of the continuous system, the rate of 
interactions per picosecond increases to around 4.7 x 10 4 
coll/ps and the relative efficiency of the DMD simulations 
drops to 1.25. Of course these results depend strongly on 
a number of conditions, including the physical conditions 
of simulation, and are particularly sensitive to the den- 
sity. For example, for a gaseous system where T* = 10 
and p* = 0.01, the DMD simulation with AV = kT/2 
becomes more than 120 times more efficient than the 
SMD simulations. From an algorithmic point of view, 
the low density DMD simulations resemble an adaptive 
time step simulation approach in which large time steps 
are utilized to integrate the equation of motion of a given 
particle when it is not interacting with any others, and 
when another particle is encountered, the time step is 
reduced to properly evolve the system through the inter- 
action region. The energy-terracing scheme and event 
detection algorithm naturally provide an adaptive ap- 
proach in which an interacting pair may experience many 
events in rapid succession as the pair evolves through an 
interaction region. 

VI. CONCLUSIONS 

In this article a general framework for performing 
event-driven dynamics of rigid bodies has been presented 



that makes use of a mapping of a continuous potential 
onto discrete values. The evolution of the discretized- 
energy system consists of periods of free propagation of 
the system punctuated by impulses generated at discrete 
times that correspond to moments when the underlying 
continuous interaction potential for a pair of particles 
hits a critical value. An adaptive grid method using cu- 
bic interpolation was presented to facilitate the search for 
the earliest interaction time. The effect of the impulses 
on the subsequent dynamics of the system was analyzed 
using conservation conditions, resulting in a simulation 
algorithm that solves the evolution of the system within 
numerical precision. The trajectories are time-reversible, 
and exactly obey all applicable conservation laws by con- 
struction. 

The method was demonstrated on a system of dense 
rigid rods interacting by a discretized version of a sim- 
ple modified Lennard-Jones pair interaction potential in 
which the effective distance depends explicitly on the rel- 
ative orientation of the rods. Although static correla- 
tions were reasonably well represented in the discretized 
potential system provided the energy mapping was not 
too coarse, dynamical correlations and related quantities 
such as the orientational relaxation time depend sensi- 
tively on the level of discretization of the continuous po- 
tential. It was found that potential energy steps on the 
order of kT/2 were required to reproduce results from 
simulations of the continuous potential system. 

It should be emphasized that the implementation 
tested here, in which an evenly-spaced energy discretiza- 
tion was used, is the simplest choice of mapping of the 
continuous system onto a set of discrete energy levels. It 
is quite likely that some other distribution of energy levels 
would lead to substantial improvements in both the qual- 
ity of the simulation results and the overall efficiency of 
the algorithm. However such modifications to the imple- 
mentation must be carefully considered, since typically 
one is only interested in rough qualitative behavior of a 
system that is not overly sensitive to details of an inter- 
action potential. In this case, the goal is to construct a 
simple model that demonstrates the relevant physics. A 
few, well-chosen potential energy discontinuities can be 
expected to meet this goal. 
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